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Abstract. We consider the three-dimensional (3D) mean-field model for the Bosc- 
^ ' Einstein condensate (BEC), with a ID nonlinear lattice (NL), which periodically 

changes the sign of the nonlincarity along the axial direction, and the harmonic- 
oscillator trapping potential applied in the transverse plane. The lattice can be created 
as an optical or magnetic one, by means of available experimental techniques. The 
I ■ objective is to identify stable 3D solitons supported by the setting. Two methods 

are developed for this purpose: The variational approximation, formulated in the 
framework of the 3D Gross-Pitaevskii equation, and the ID nonpolynomial Schrodingcr 
O . equation (NPSE) in the axial direction, which allows one to predict the collapse in the 

framework of the ID description. Results are summarized in the form of a stability 
region for the solitons in the plane of the NL strength and wavenumbcr. Both methods 
^ , produce a similar form of the stability region. Unlike their counterparts supported by 

the NL in the ID model with the cubic nonlinearity, kicked solitons of the NPSE 
■ cannot be set in motion, but the kick may help to stabilize them against the collapse, 

by causing the solitons to shed excess norm. A dynamical effect specific to the NL 
is found in the form of freely propagating small-amplitude wave packets emitted by 
. perturbed solitons. 



X 



Bose-Einstein condensates in nonlinear lattices 



2 



1. Introduction and the model 

The use of periodic potentials, induced by optical lattices, for steering matter waves in 
Bose-Einstein condensates (BECs) is a vast research cirGci, £is demonstrated in reviews p]- 
[4]. An important aspect of this topic is that the lattice potentials, balancing the cubic 
nonlinearities induced by inter-atomic collisions in the BEC, help to create and stabilize 
solitons. In particular, the lattices play a critical role in stabilizing two-dimensional 
(2D) solitons against the collapse in the condensate with intrinsic self-attraction [5]. 

Theoretical and experimental studies of the soliton dynamics in periodic potentials 
were recently extended for nonlinear lattices (NLs), which may be induced by a 
spatially periodic modulation of the local strength of the nonlinearity (the respective 
effective nonlinear potentials are often called pseudopotentials [6]). In BEC, the spatial 
modulation can be implemented by means of the Feshbach resonance controlled by 
properly patterned external fields (in the quasi-lD BEC, a combination of linear and 
nonlinear lattices may also be induced by periodically modulating the strength of 
the tight transverse confinement in the axial direction [7]). As well as their linear 
counterparts, NLs have drawn much attention in connection to their potential for the 
creation and control of matter- wave solitons in a number of different settings, see original 
works [8]- [18] and review [I]. While NLs readily support stable ID solitons [Ej-[T7], it 
has been found difficult, albeit sometimes possible, to stabilize 2D solitons against the 
collapse by means of the NL-induced pseudopotentials [18]. As concerns the ID solitons 
in models with the cubic nonlinearity, the numerical analysis also reveals that they 
feature mobility in the presence of NLs [H [9] . 

Thus far, no example of stabilization of 3D solitons by NLs has been reported. 
On the other hand, the creation of multidimensional solitons can be facilitated by 
combinations of linear and nonlinear lattices, as shown, in particular, in Ref. [TT] . 
where 2D solitons supported by crossed ID lattices, one linear and one nonlinear, were 
reported. Still earlier, it was demonstrated that the ID linear-lattice potential, acting 
together with periodic temporal modulations of the nonlinearity, which may be induced 
by the Feshbach resonance controlled by a time-periodic external field, can stabilize 3D 
solitons [19]. 

A natural possibility for the creation of 3D solitons is suggested by a combination 
of the ID NL with the usual harmonic-oscillator linear trapping potential acting in 
the plane perpendicular to the NL axis (similar settings, but with linear ID lattices, 
were shown to be efficient in the creation of 3D gap solitons in the BEC with the 
repulsive intrinsic nonlinearity [20]). This setting is the subject of the present paper. 
We tackle it by means of two different methods, namely, the variational approximation 
(VA) applied to the underlying 3D Gross-Pitaevskii equation, and, in a more accurate 
form, the effective ID nonpolynomial Schrodinger equation (NPSE), which was efficient 
in description of many other settings dominated by the interplay of the tight 2D 
confinement and nonlinearity [21], [22], [7], including the action of linear lattice 
potentials in the axial direction [23]. 
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Thus, we consider a dilute BEC of bosons with mass m confined in the transverse 
plane by the isotropic harmonic-oscillator potential with frequency cjj_, V(x, y) = 
(l/2)moj 2 _ (x 2 + y 2 ). The corresponding 3D Gross-Pitaevskii equation is rescaled by 
measuring time and coordinates in units of uj 1 and the transverse-confinement radius, 
a± = \Jhj (mwi), respectively (hence the energy is measured in units of Huj±): 



dijj 



- l -V 2 + l -{x 2 + y 2 )+2ng(z 



4>, (1) 



with the condensate's wave function normalized to unity, 

\^{r,t)\ 2 al 3 r = 1. (2) 

The interaction strength in Eq. (JTJ is 

g(z) = 2(N-l)a s (z)/a ±J (3) 

where N is the number of atoms and a s the ^-dependent s-wave scattering length of the 
inter-atomic potential, the NL corresponding to the periodic dependence, 

g(z) = g + gi cos (2k z). (4) 

Below, we consider the most fundamental case of go = 0. Placing the center of the 
soliton at z = 0, we assume g± < 0, to support the soliton by the locally attractive 
nonlinearity. 

The ID periodic modulation of the local scattering length, implied by Eqs. (jlj) and 
([3]), can be implemented in an optical lattice, produced by the interference of a pair of 
counterpropagating laser beams controlling a s via the Feshbach resonance [21]. In that 
case, the period of the resulting NL is limited by diffraction to > 1 /im. More often, the 
Feshbach resonance in experiments with BEC is controlled by the magnetic field [25J. 
In that case, the ID periodic structure can be built as a magnetic lattice, imposed by a 
properly designed set of ferromagnet films [26], with the respective fabrication limit on 
the NL period also amounting to > 1 /im. Then, assuming that the trapping potential 
confines the transverse size of the condensate, as usual, to the same order of magnitude 
(a few microns), one may conclude that the solitons are built of several thousand atoms 



The rest of the paper is organized as follows. The direct VA is developed in Section 
II, and the effective ID NPSE is derived in Section III. In both cases, we find the 
stability region for solitons in the plane of (k, \gi\). The mobility of the solitons is 
tested in Section III by means of direct simulations of the evolution of axially kicked 
solitons, in the framework of the NPSE (it is concluded that the solitons are not mobile 
in the present setting; however, the kick may effectively stabilize solitons against the 
collapse). The paper is concluded by Section IV. 
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Our first aim is to apply the VA to Eq. (ED), following the lines of Ref. [31]. To this end, 
we notice that Eq. ([!]) can be derived from the Lagrangian density, 



C 



i dip dip*. 



_|VVf --(* 2 + y 2 )H 2 -7r^)H 4 , 
and make use of a time-dependent Gaussian ansatz, 



exp 



Pj(r,t) 







z 2 




a'i(t) 


°f(*) 



(5) 



(6) 



where r\ = x 2 + y 2 , and cr\\(t) and f3±(t), (3\\(t) are time- dependent variational 

parameters. This wave function is an exact one for non-interacting bosons (g = 0) in 
the harmonic trap. 

Inserting the ansatz into Lagrangian density ([5]) and performing the spatial 
integration, we arrive at the effective Lagrangian, 



+ 0\\ 



1 

2 L 
o-l 



1 



+ \o\$i + c-1) 



2a 



£0 + 9ie 



-k 2 a 2 /2 



(7) 
(8) 



v2vr o-\a\\ 

with the overdot standing for time derivatives. The respective Euler-Lagrange equations 
take the form of 



fil 



C7_l + a± 



2a± ' 

_ 3_ 

2cr„ ' 



-k 2 a 2 /2 



27T cr^cj|| 



9o + 0ie 



-fcV 2 /2, 



1 + k 2 a 2 ' 



(9) 
(10) 

(11) 
(12) 



Equations Qj and ffTUj) show that, as usual, chirps and are determined by the 
time dependence of a± and cry, while Eqs. (TTT]) and ( |T2l) correspond to the equations of 
motion of a mechanical system with two degrees of freedom, whose energy is 



E 



1.2 , 1.3 



2 

U(a±,a\ 



1 2 1 

2^ + 2aJ 
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go + gie 



-k 2 <T 2 /2 



4a 



(13) 



2V27T a]_a\\ 

Next, we look for stationary configurations corresponding to minima of potential 
energy ( fT3"j) . by demanding dU/da± = dU/da\\ = 0, which yields 

. ~k 2 a 2 /2 

i . go + gie »' ( . 

o-± = — + 7=-^ , (14) 
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=-3 + 7^r-o -■ (15) 



<7jj V27r cr^ajj 

These equations, which are fixed points of Eqs. ( TTTT) and ( TT2|) . with <t_|_ = cry = 0, 
can be solved numerically. The solutions provide for a minimum of the energy under 
the necessary condition that the Gaussian curvature Kq of energy surface U(a±,a\\) is 
positive, i.e., 

* ( *£-V > o . (16) 



da]_ da 2 \ da± da u r 

Further, low-energy excitations of the condensate around the stationary state are 
represented by small oscillations of variables o~j_ and cry around the equilibrium point 
defined by Eqs. (]14p and ffl5|) . The calculation of the corresponding normal- mode 
frequencies, Q, is thus reduced to finding eigenvalues of the respective Hessian matrix, 

/ d 2 U/da\ d 2 U/da ± daA 

y d 2 U/da ll da ± d 2 U/da 2 J ' K '> 

while the associated mass matrix is 



M 





l 

2 



Then, the eigenfrequencies are found as the solutions of equation 

det (A - tt 2 M) = . (18) 

Widths <Tj_ and a» of stable solitons, as found from the numerical solution of 
Eqs. and f fT5l) . are plotted in the left panels of Fig. [1] as functions of the NL 

strength parameter, [recall we set go = and g\ < in Eq. (TJJ], along with 
the widths obtained by solving numerically the ID NPSE (see below). The comparison 
demonstrates close proximity of the predictions of the VA to the results produced by the 
NPSE. In the right panels of Fig. [TJ we plot eigenfrequencies Q\ and Sl 2 of excitations 
around the bright soliton, as found from Eq. (fTSl) . 

In Fig. H] we plot the stability diagram in the plane of (k, \gx\) for the bright 
solitons trapped in the NL with wavenumber k and strength \gi\ (again, with g = 0). 
The stability region, defined as that where solutions of Eqs. ( fT4l) and ( fl5l) yield energy 
minima, is bounded by the dashed lines (strictly speaking, the bright solitons are only 
metastable, because for g = and g\ < the true ground state is the collapsed one, 
with potential energy U = — oo). Below the lower dashed line, the soliton is subject 
to spreading along longitudinal axis z (which may be considered as a manifestation of 
the derealization transition, which was earlier studied in linear lattices [30]). Above 
the upper dashed line, the soliton is destroyed by the collapse. Very close to the upper 
dashed line, the VA predicts a bistability, i.e., coexistence of two stable solitons at the 
same values of k and \g\\ (in Fig. dj only the solitons with the lowest energy are shown 
in the case of the bistability). In Fig. [2] we also plot the respective stability region 
(between the solid lines) as obtained from the numerical solution of the NPSE (see 
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Figure 1. (Color online) Left panels: transverse and axial widths, <r± and a\\ (dashed 
and solid lines, respectively) of stable bright solitons, as functions of the NL strength, 
\gi\, with go = and g\ < 0, obtained from the variational approximation. Filled 
circles and squares depict the widths obtained from a numerical solution of the NPSE. 
Right panels: two eigenfrequencies f2i and (dashed and solid lines) of collective 
excitations around the stable bright soliton, vs. \gi\, as obtained from the variational 
approximation. The results are shown for three values of the NL wavenumber: k = 0.2 
(upper panels), k = 0.5 (middle panels), and k = 0.8 (lower panels). 



2 - 




Figure 2. (Color online) The stability diagram for the solitons in the plane of 
wavenumber k and strength \gi\ of the NL (with go = 0). The solitons are stable 
between the dashed lines, according to the variational approximation, and between the 
solid lines, according to the NPSE. The dot-dashed line is the lower bound predicted 
by the one-dimensional cubic Gross-Pitacvskii equation. 

below). The shrinkage and disappearance of the stability region at large values of k is 
quite natural, as in that case the rapidly oscillating nonlinearity in Eq. (pjj tends to 
average itself to zero. 

The prediction of the usual ID cubic Gross-Pitaevskii equation is also plotted in 
Fig. |2j In that case, there is only one stability boundary (the dot-dashed line), as the 
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ID equation with the cubic nonlinearity does not predict the collapse. Accordingly, the 
above argument concerning the disappearance of the stability region at large k does 
not apply to the cubic equation, because, in the limit of the small NL period, %/k, 
the soliton may compress itself into a single potential well, and this trend will not be 
aborted by the onset of the collapse. 



3. The nonpolynomial Schrodinger equation (NPSE) 

3.1. The derivation and imaginary-time evolution 

A more accurate description of the solitons is provided by the NPSE, which can be 
derived by means of the semi- variational approach from the full 3D equation (JTJ), using 
the method developed in Ref. [21]. To this end, we adopt the following ansatz, which, 
unlike the above one ([6]), contains arbitrary functions of the longitudinal coordinate, 
o~(z,t) and f(z,t), accounting for the transverse width and axial wave function of the 
condensate: 

1 



7jj(r,t) 



exp 



y/ira(z,t) 

Note that, as follows from Eqs. (fT9"]) and 
1: 



/CM). 



(19) 



N- 



1D 



\f{z)\ 2 dz = l. 



x 2 + y 2 

2 KM)) 2 . 

(J2J) , the norm of the axial wave function is also 

(20) 



Substituting ansatz (lT9~j) into Lagrangian density fl5}, performing the integration over x 
and y, and omitting spatial derivatives of the transverse width (this corresponds to the 
adiabatic approximation, which is known to produce accurate results in other settings 
[2~T]). we derive the respective Lagrangian density, 

2 1 / 1 1 I fl 4 

l/l 2 -^)^- (21) 
2 cr z 



C 



2 U dt 



f 



dp 



dt 



df 



dz 



G- 



Varying it with respect to f*(z, t) and a(z, t) gives rise to the system of Euler-Lagrange 
equations: 



.df 



a 



dt 

1 



l/p 



Id 2 1/1 9 \ 

1 + <?(2)l/| 2 , (23) 

Inserting Eq. ( 123]) into Eq. ( 122]) . we obtain a closed- form equation for the axial 
wave function, which is tantamount to the NPSE derived in Ref. [21], but with the 
^-dependent interaction strength, g(z): 



.df 



2dz* 



+ (3/2)g(zM 
V± + 9(z)\f\ 2 



f ■ 



(24) 
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In the weak-coupling regime, i.e., \g(z)\ \f(z, t)\ 2 <C 1, one can expand the 
nonpolynomial term in Eq. ( jMJ) , arriving at the cubic-quintic nonlinear Schrodinger 
(with term 1 representing here the transverse ground-state energy), 



.df 



^ + i + ^)l/l 2 + ^) 2 l/l 4 



The cubic-quintic nonlinearity for the tightly confined BEC was also derived by means 
of different approaches [29]. On the other hand, in the strong-coupling regime, 
g(z) \f(z, t)\ 2 ^> 1 (which is relevant only for the repulsive sign of the nonlinearity, 
g > 0), the NPSE amounts to the nonlinear Schrodinger equation with the quadratic 
nonlinearity (see, e.g., Ref. [32]): 

.df 

l Tt = 



Here we aim to analyze bright solitons within the framework of Eq. ()24p with 
g(z) given by Eq. (jlj), with g = and g\ < 0, as said above. Results were obtained 
from numerical solutions based the finite-difference Crank-Nicolson predictor-corrector 
algorithm [33] . 

First, by simulating the NPSE in imaginary time, we study the formation of bright 
solitons. In particular, at g\ = —0.4, the soliton does not self-trap, slowly degenerating 
towards a uniform configuration along axial direction z. Instead, at g\ — — 1 the bright 
soliton self-traps quickly, representing the ground state of Eq. ( 124]) . In Fig. [3] we 
plot typical examples of the axial density, p(z), of the so obtained stable bright soliton 
trapped in the NL, comparing the NPSE results to those predicted by Gaussian ansatz 
OH]) (the stability of the solitons was verified by real-time simulations, see below). The 
figure shows that the stable solitons are localized around one potential minimum of the 
NL (at x = 0). The NPSE profiles are quite close to their variational counterparts, see 
also the left panels in Fig. [TJ The main quantitative, although not very large, difference 
between the VA and NPSE is the prediction of the critical strength for the onset of the 
collapse, as seen in Fig. |2] 



4. Real-time dynamics 

4-1. Stability of the solitons 

The next step is the study of the real-time dynamics of the quasi-lD BEC trapped in the 
NL. First of all, simulating NPSE ( 124]) in real time, we have found that all the existing 
solitons are stable, with initially perturbed wave functions featuring small oscillations 
around the solitonic configurations. As shown in Fig. H] an interesting dynamical feature 
is observed if the initial wave function is the Gaussian with a width close to that of the 
soliton: expulsion of two small waves from the Gaussian peak (which is represented by 
the column centered at z = in Fig. H]). The emitted waves rapidly move in opposite 
directions, while the remaining central peak relaxes into a stationary soliton. This effect 
is interesting because the same is not observed in linear lattices, which would readily 




Figure 3. (Color online) Typical examples of axial density p{z) = |/(z)| of stable 
solitons. The solid and dashed lines display the results produced by the NPSE and 
variational approximation based on ansatz ([6]), respectively, for three different values 
of the interaction strength |<7i|, fixing go = and k = 0.5. Here and in Figs. IH[5j and 
|6l the green sinusoidal line represents the periodic modulation function of the local 
nonlinearity defined in Eq. (|4]). 
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Figure 4. (Color online) Dynamics of a Gaussian wave packet with the initial shape 
close to the ground-state soliton. The axial density, p{z), is displayed at different 
values of real time t, as obtained from simulations of Eq. ([24]) . At t = (the initial 
condition, not shown here), there is only the Gaussian centered at z = 0, with axial 
width 1.7. The parameters are go = 0, g\ = —1.2, and k = 0.5. 



trap the radiation "garbage" emitted by the central peak, while the NL is not felt by 
the small-amplitude waves, hence they may escape freely. 

4-2. Immobility of the trapped solitons 

The mobility of solitons trapped in the NL can be tested by applying a kick to 
initially quiescent solitons [8]. For this purpose, Eq. ( 124")) was simulated with initial 
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Figure 5. (Color online) The evolution of the kicked soliton with initial velocity 
v = 0.4. Axial density p(z) is plotted at different values of real time t, as obtained 
from simulations of Eq. (|24|) . Here, go = 0, <?i = —1.2, and k = 0.5 are fixed. 



condition f(z, t = 0) = f so i(z) e roz , where v is the magnitude of kick, i.e., the initial 
velocity imparted to soliton f so \(z), which was produced by means of the imaginary-time 
simulations of the same NPSE. To present a typical result, we fix g — 0, g% — —1.2, 
and k = 0.5, and perform real-time simulations at increasing values of v. As shown in 
Fig. [51 at v — 0.4 we observe ejection of small-amplitude waves from the soliton (cf. 
Fig. H|), while the central peak remains trapped at the initial position, relaxing back 
into a stationary soliton, with a somewhat smaller value of the norm. 

With the increase of v, the amplitude and the velocity of the ejected waves increases, 
but the remaining soliton stays put. This is a noteworthy difference from the soliton 
dynamics in ID NLs with the cubic nonlinearity, where the soliton may be set in motion 
by the kick [8]. Eventually, if the kick is too strong, it destroys the soliton. In particular, 
for \gi | = 1.2 the destruction is observed at v > 0.48, see an example in Fig. [6]for v = 0.6. 

The critical velocity, v c , at which the kicked soliton is destroyed, is shown in Fig. [7] 
as a function of the NL strength, \gi\. The figure features a linear growth of v c with NL 
strength \gi\, at sufficiently large values of \g\\. This fact can be explained by estimating 
the critical velocity as that at which the respective kinetic energy of the kicked soliton, 
(l/2)M so \v 2 , is equal to height Vpn of the effective Peierls-Nabarro potential induced 
by the nonlinear (pseudo) potential. The mass of the soliton, M so \, is proportional 
to its norm, which is 1, according to Eq. (TSUI) . Further, the potential-energy density 
corresponding to Eq. ( J24|) actually coincides with the potential part of Lagrangian 
density (I2T1) . that should be evaluated with the help of expression (|23|) for a. Then, 
a straightforward consideration of Eqs. (|2T|) . (123]) and (1241) yields the following scaling 
relations in the limit of large \g\\: |/(0)| 2 ~ a^ 1 ~ W~ l ~ |<7i|~ , where W is the 

2 m 

axial size of the soliton, and, eventually, Vpn ~ |gi| • Thus, the threshold condition, 
(l/2)M so \v 2 = Vpn, explains the linear proportionality between v c and \gi\, which is 
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Figure 6. (Color online) The same as in Fig. [5j but for the soliton kicked with initial 
velocity v = 0.6. 
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Figure 7. The critical velocity, v c , for the destruction of the kicked soliton versus the 
strength of the nonlinear lattice, \gi\. Here, go = and k = 0.5 are fixed, as before. 

observed in Fig. at large \gi\. The vanishing of v c at \gi\ m 0.75 in Fig. [7] is a 
consequence of the fact that the soliton with this value of \g±\ lies at the edge of the 
triangular area in Fig. El i.e., it does not exist as a stable mode even without being 
kicked. 

The fact that the kick induces emission of radiation from the soliton may be used 
to stabilize them in the collapse domain of Fig. [2j To this end, we take, for example, 
the Gaussian wave packet with g± = —1.5 and k = 0.5, which falls into the region of the 
collapse. Real-time simulations of Eq. ( 124)) demonstrate that, if the Gaussian is kicked 
hard enough, it does not blow up, but rather forms a stable soliton, in a combination 
with the emission of radiation waves. This is shown in Fig. [8] for initial velocity v = 0.4. 
The initial configuration evades the blowup because the emission of radiation reduces 
the norm of the remaining soliton, pushing it beneath the collapse threshold, see the 
lower panel in Fig. [SJ This observation suggests that the relaxation of the perturbed 
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Figure 8. (Color online) Upper panel: The dynamics of the Gaussian wave packet, 
kicked with initial velocity v = 0.4, whose parameters belong to the collapse domain in 
terms of Fig. [2] gx = —1.5, k = 0.5, and go = 0. At t = (the initial condition shown 
by the solid line), there is only the column centered at z = 0, which represents the 
Gaussian of axial width 1.7. Lower panel: the norm of the wave function in interval 
— 15 < x < 15, as a function of time. 

soliton via the emission of radiation proceeds faster than the onset of the collapse, which 
attests to the robustness of the solitons. Finally, as expected, if the kick is too hard (in 
the present case, this means v > 0.6), it destroys the Gaussian wave packet, causing its 
complete decay into radiation. 

5. Conclusions 

We have reported results for 3D matter-wave solitons supported by a combination 
of the axial ID NL (nonlinear lattice), which periodically reverses the sign of the 
nonlinear interaction, and the tightly trapping harmonic-oscillator potential acting in 
the transverse plane. The results were obtained by means of two distinct approaches: 
The VA (variational approximation), which was applied directly to the 3D Gross- 
Pitaevskii equation, and the ID NPSE, that was derived from the underlying 3D 
equation. Previous works did not study the stabilization of 3D solitons by NLs. The 
main result, produced by means of both methods in similar forms, is the stability domain 
for solitons in the plane of the NL strength and wavenumber. The usual ID cubic Gross- 
Pitaevskii equation with the NL cannot produce adequate results, as it does not give 
rise to the collapse, which is the most important stability-limiting factor. 

Another essential difference of the solitons produced by the NPSE with the NL 
from their counterparts in the case of the cubic NL is that the solitons are immobile 
in the framework of the NPSE: The kick applied to the soliton either leaves it pinned, 
or, eventually, destroys it. The critical size of the kick which destroys the soliton was 
found to be proportional to the strength of the NL, provided that the strength is large 
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enough; an explanation for this dependence was proposed. On the other hand, the kick, 
if applied to the wave packet created above the collapse threshold, may help it to shed 
off the excess norm and thus stabilize itself against the collapse. A related dynamical 
effect, which demonstrates the difference of the NLs from linear lattices, is that wave 
packets relaxing into solitons can emit small-amplitude waves, which freely propagate 
in the system. 

A challenging extension of the analysis may be to develop it for the setting with a 
2D NL and the ID trapping potential acting in the transverse direction. The 2D version 
of the NPSE was developed previously, but in the absence of the NL [2TJ EI]. In this 
case, one may expect the existence of both fundamental and vortical 2D solitons. 
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